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ABSTRACT 

It  Is  suggested  here  that  the  electronic  wave  functions  of  one 
electron  in  the  field  of  two  nuclei  may,  in  some  cases,  provide  a  good 
basic  set  of  one-electron  orbitals  for  variational  calculations  of 
electronic  wave  functions  of  several-electron  diatomic  molecules. 
General  formulas  and  computational  methods  for  the  necessary  electronic 
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integrals  appropriate  to  the  homopolar  diatomic  molecule  is  described. 
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NYO-10,423 

ON  THE  USE  OF  DIATOMIC  ORBITALS  IN  CALCULATIONS 
OP  THE  ELECTRONIC  WAVE  FUNCTIONS 
OF  DIATOMIC  MOLECULES 

James  W.  Cooley 

1.    Introduction. 

One  approach  to  the  problem  of  calculating  the  electronic 
wave  functions  of  molecules  is  to  employ  variational  methods 
with  N-electron  trial  wave  functions  constructed  from  a  basic 
set  of  one-electron  functions  or  orbitals.   The  Slater-type 
atomic  orbitals  (STAO's), 

where   r,  0,  0  are  polar  spherical  coordinates  and  the  Y.  , 

are  spherical  harmonics  centered  on  each  of  the  nuclei,  have 

2 
been  used  quite  extensively.    More  recently,  in  the  case  of 

diatomic  molecules,  orbitals  of  the  form 

1_ 
^PS  "  [(A^-1)(1-pl^)  ]^    x\^exp{-a}^-^[l+lm0)  (1.2) 

where   A,  |x,   and  0     are  prolate  spheroidal  coordinates  with 

3  if  c,  6 
the  nuclei  as  foci,  have  been  found  '  '"•*    to  present  certain 

advantages  in  the  calculation  of  the  required  energy  integrals 

and  in  the  accuracy  of  the  results  obtained. 

It  is  suggested  here,  that  in  some  cases  it  may  be  advantageous 

to  use,  instead,  diatomic  orbitals   (DO's),  i.e.  the  solutions  of 

one-electron  diatomic  molecules  as  the  basic  orbitals. 

7 
In  1951  Hylleraas   calculated  an  approximate  solution  for 
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H   by  assigning  the  two  electrons  to  the   ISc   solutions  of 
the  H  "^  -  like  molecule  with  nuclear  charges  of   1   for  one 
electron  and  l/2   for  the  other.   He  estimated  the  electronic 
repulsion  energy,  however,  by  using  only  the  first  terms  of  the 
two  Infinite  expansions  occurring  in  the  solutions.   Shortly 

o  , 

after,  Recknagel   used  H^  -  like  solutions  for  the  outer 

9 
electrons  of  N  .   Then,  much  later,  in  195^,  Wallls  and  Hulbert 

extended  the  Hylleraas  calculations  by  including  the  first  two 

terms  in  each  of  the  above-mentioned  expansions.   They  also 

discussed  some  of  the  problems  Involved  in  using  the   DO=s 

for  molecules  larger  than  H^ .   Rahman   used  orbitals  approximately 

of  the  form  of  the  H  "*"  solutions  but  his  treatment  was  actually 

only  a  limited  application  of  the  basic  set   (l,2). 

The  wave  equation  for  the  DO's   is  separable  in  prolate 

spheroidal  coordinates  and  the  exact  DO's   in  these  coordinates 

may  be  expressed  in  the  form 

-  m 
Xj^Q  =    [(A^-l){l-H^)]^   g(A)f(M.)  exp(-pA-?M.+im0)       (l,3) 

where  p   is  related  to  the  electronic  energy  and  p  =  +  p  or   0, 
depending  on  the  choice  of  several  alternatives  for  representing 
the  solution  of  the  separated  |x  equation.   The   p  =  0  form 
can  arise  only  for  the  homopolar  molecule.   The  functions   g(A) 
and  f(iJ.)   each  may  be  represented  by  infinite  expansions  in 
functions  of   A  and  \i,      respectively,  whose  coefficients  can 
be  calculated  by  methods  to  be  described  below.   An  immediate 
consideration  which  comes  to  mind  when  considering  the  ^j^q's 
as  an  alternative  to  the  X  Q-p      is  how  close  the   ^do's   come 
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to  lying  in  the  linear  manifold  spanned  by  the  4op,    which 
one  would  use  in  a  mulecular  calculation.   This  would  be  equivalent 
to  being  able  to  represent  g(A)   and  f(iJ,)   as  polynomials  in   A 
and  \i,      respectively.   On  examining  the  differential  equations  for 
g(A)   and  fda),   it  is  found  that,  in  the   p  =  0  case,  f{\x) 
admits  a  power  series  expansion  in  ^i  but  gCx)   does  not*.   This 
implies  that  one  cannot  approximate  g(A)   uniformly  by  polynomials 
in   A   of  Increasing  degree.   Furthermore,  one  cannot  say  that  any 
molecular  wave  function  obtained  with  '^t-jq,    can  be  surpassed  in 
accuracy  by  taking  sufficiently  many   Aqp   iri  the  basis  set. 
This  argument  does  not  support  a  preference  for   A^q,    but  does 
show  that  they  are  sufficiently  different  from  Zor^,         to  warrant 

or  S 

a  computational  investigation  of  their  usefulness  in  molecular 
calculations . 

One  of  the  more  immediate  objectives  of  the  planned  calculations 
with  DO's  will  be  to  see  if  a  given  degree  of  accuracy  can  be 
attained  with  a  smaller  basic  set  of  orbitals  when  using   Anoi 
instead  of   Aaoic;   °^  /^PS'<^"   '^^is  will  be  quite  important 
in  configuration  interaction  calculations  since  the  number   of 
configurations  increases  extremely  rapidly  with  the  number  of 
orbitals  used.   Harris   has  pointed  out  a  number  of  advantages 
in  using  the  spheroidal  coordinate  orbitals  rather  than  STAO's 
which  also  apply  to  DO's.   One  of  the  most  important  is  that 
with  a  single  generic  form  for  all  orbitals,  rather  than  having 
orbitals  on  the  different  centers  of  the  molecule,  a  more  systematic 
computational  procedure  may  be  used  for  all  integrals.   One  further 

*  See  the  discussion  on  page  620  of  reference  11. 
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advantage  Is  that  orbltals  In  prolate  spherical  coordinates  are 
not  subject  to  the  overcompleteness  caused  by  using  large  numbers 
of  STAG'S  on  several  centers  of  a  molecule*.   It  will  be  of 
particular  interest  to  compare  predicted  local  properties  of  a 
molecule.   For  example,  experience  has  shown**  that   STAG  wave 
functions  give  too  low  a  value  of  the  electric  field  gradient 
at  the  nuclei  as  compared  with  wave  functions  in  prolate 
spheroidal  coordinates . 

In  the  next  two  sections,  the  discussion  and  the  formulas 
derived  are  general  enough  to  apply  to  the  heteropolar  molecule 
and  to  any  representation  one  chooses  to  use  for  the  separated 
solutions  of  the  one-electron  wave  equation  in  prolate  spheroidal 
coordinates .   The  machine  program  referred  to  and  the  computations 
performed  for  the  homopolar  molecule  are  somewhat  more  specialized 
and  are,  therefore,  described  in  the  appendices. 

2.    The  Diatomic  Orbltals 


The  wave  equation,  in  atomic  units,  for  an  electron  in  the 
vicinity  of  two  fixed  nuclei  A  and  B,   at  a  distance  R  apart 

and  having  nuclear  charges   ^   and  11,,  is 

3,  D 


(^  A  +  C  r  "^  +  Zy,r-^   +   b)X  =  0'  ^2.1) 


'a  a     ^b  b 

where  A   is  the  Laplacian  operator,   r   and  r,   are  the 

*  It  is  possible  to  encounter  overcompleteness  in  a  set  of 
spheroidal  coordinate  functions  for  some  values  of  the  scale 
parameters  when  these  are  varied  independently. 

**See  page  150^  reference  12.   Also,  see  reference  15- 
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distances  of  the  electron  from  nuclei  A  and  B,   respectively, 
and  the  eigenvalue   e   is  the  energy  of  the  electron.   If 
expressed  in  the  prolate  spheroidal  coordinates 

A  =  (r^  +  r^)R"\ 

1  (2.2) 

and  ^,   the  aximithal  angle  of  the  electron  about  the  molecular 
axis,  (2.1)  can  be  separated  and  its  solution  can  be  written 

X  =  M(M.)A(A)exp(imj2f)  (2.3) 

where  m  =  0,  +  1,  +  2,  ...   and  where  M{|j.)   and  /\(a)   are 
solutions  of  the  differential  equations 

[£^(n)  -  p^n^  +  2p^n  +  A]M(m.)  =  0  (2.4) 

[j:^(a)  -  p^A^  +  2P2A  +  a]A(a)  -  0.  (2.5) 

A  is  the  separation  constant  and 

Pi  =1  ^^^a  -  ^b^^  (2.6) 

P2  -I  RUa  +  Cb).  (2.7) 

p2  =  -  I  R^E.  (2.8) 

The  differential  operator   <C  (z)   is  defined  by 

£  (z)  =  ^  (z^  -  1)  ^^ 2_  (2.9) 

^  dz  dz   z  -1 

and  has  the  eigenvalues  £{£+!)      i  =  m,  m  +  1,  ...  .   The 
corresponding  eigenfunctions  of  ^^^^"^      ^^^  ''^he  associated 
Legendre  functions   ^^£(2)   and  q|^(z).   Without  loss  of  generality. 
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we  restrict  ourselves  to  positive  values  of  m  for  the  present. 

Equations  (2.4)  and  (2.5)  are  of  the  same  form  with  the  only 
difference  being  In  the  ranges  of  the  Independent  variables; 

-  1  ^  [X  <  1 

(2.10) 

1  ^  A  <  00  . 

A.  H.  Wilson   and  Baber  and  Hasse   obtained  and  Investigated 

various  forms  of  the  solutions  of  (2.4)  among  which  one  finds 

00     „ 
M(n)  ^-^  F^vliyi)  (2.11) 


k  k' 


k=m 


TT  m   .      QD 


M(m,)  =  (1-^2)2    e±   P^  Y^   f^(l^)^  (2.12) 

k=0 

where  the  F;,'s   of  (2.II)  satisfy  a  three-term  recurrence  relation 

If  p   =  0,   I.e.  for  the  homopolar  molecule,  and  the   ^j^'s   of 

(2.12)  satisfy  a  three -term  recurrence  relation  for  arbitrary 

--14 
p, ,  I.e.,  for  the  heteropolar  molecule.   Baber  and  Hasse   point 

out  that  (2.12)  converges  rapidly  for  large   p  and  Is  suitable 

for  the  entire  range  of  p.   However,  for  moderate  p  and  for 

P-,  =  0,   (2.11)   Is  more  convenient  to  use. 

Wilson  did  not  obtain  a  satisfactory  solution  of  (-2.5)  and, 

before  questions  of  the  existence  and  nature  of  such  solutions 

1^  7 

were  adequately  treated  -^ ,   Hylleraas   used  the  form 

(A)  =  (a2  -  1)2  ""  e-P^^-1^  ^  E^   L^^"^^(2p(A-l);, 

k=0 

(2.13) 

where  the   LJ''^  's   are  the  associated  Laguerre  functions. 
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Jaffe   claimed  erroneously  that  the  series  (2.13)  did  not 
converge*   and  obtained  a  solution  In  the  form 

1 
A(A)  =  (?2_i)2  ""  e-P^(A+lf  ^  G^Z^  (2.14) 

k=0 
where 

C  =  p^P    -  m  -  1 

Z  =  (A-1)/(A+1) 

In  both  (2.15)  and  (2.l4)  three-term  recursion  relations 
give  the  expansion  coefficients.   In  the  present  work,  Jaffe's 
solution  (2.1^)   is  chosen  only  because  it  required  less  computation 
for  the  integration  methods  used  here. 

We  will  describe  the  assignment  of  quantum  numbers  and  certain 
functional  relationships  in  a  manner  related  to  the  computational 
procedure  described  in  the  next  section.   We  first  regard  A   as 
the  eigenvalue  parameter  in  both  (2.4)  and  (2.5)  and  let   p, ,  p^ 
and  m  be  fixed.   For  each  value  of  p,   there  is  a  sequence  of 
eigenvalues   A  =  k    [p, p    ,£ ,m)      of  (2.4)  with  their  index   £ 
taking  the  values   m,  m+1,  m+2 ,  . . .   for  the  eigenfunctions 
M(ix;p,p-,  ,<g,m)   with  £-m  -   0,1,2,...    nodes,  respectively. 

Similarly,  the  eigenvalues  of  (2.5)  are  denoted  A  =  A^(p,p_,k,m) 

A     ^ 

where   k  =  0,1,2,  ...   is  the  number  of  nodes  in  the  solution 
A(A,p,Pp,k,m) .   If  we  plot   A   and  A^   on  a  graph  as  functions 
of  p,   we  obtain  a  set  of  non-crossing  A   curves  and  a  set  of 

See  reference  14,  page  570  for  a  proof  of  convergence. 
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non-crossing  A^   curves  with  each  intersection  of  an  A 
curve  with  an  A^   curve  corresponding  to  a  solution  of  (2.1). 
The  values  of  A  and  p  at  each  Intersection  may  be  written 

A  =  A{p-^,p^,n,£,m)    =   A  (p;p^,i,m)  =  A^(p;  p^,  k,m)     (2.l6) 
and 

p  =  p(p^,p2,n,^,m), 

respectively,  where 
n  =  i  +  k  +  1. 
The  corresponding  solution,  of  (2.1)  may  then  be  expressed 


X(A,ix,^;p-,,p^,n,^,m)  =  M(iJ.;p,p,  ,i,m)A(  A:p,  Po,k,m)e 


im^ 

(2.17) 
The  quantum  numbers  thus  assigned  to  the  DO's  are  the  quantum 
numbers  of  the  united  atom:  the  atom  which  the  molecule  becomes 
when  R  goes  to  zero.   It  is  Important  to  note  that  each  DO 
contains  two  parameters,   p-,   and   pp,   which  depend  on  the 
effective  nuclear  charges   ^a  and  ^b.   Like  the  orbital  exponents 
of  STAG'S,   i^a  and  ^b  may  be  regarded  as  screening  parameters 
and  treated  as  variation  parameters  in  a  variational  calculation 
of  an  N-electron  wave  function.   The  physical  interpretation  of 
DO's  is  distinguished  from  that  of  STO's,  however,  in  that  one 
now  may  think  of  the  electrons  in  a  molecule  as  occupying  orbitals 
formed  from  those  of  the  united  atom  by  having  split  the  atom  Into 
the  molecular  nuclei  rather  than  as  occupying  orbitals  obtained  by 
superimposing  the  atomic  orbitals  of  the  constituent  atoms. 
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3 •    The  Calculation  of  Electronic  Integrals 

In  this  section,  formulas  will  be  derived  for  the  two- 
electron  Coulomb  repulsion  integrals 

/;r{l);('(l)X"(2);X'"(2)r-^dT^dT2  (3.1) 

where  the  Indices   1   and  2   refer  to  the  coordinates  of  electrons 
1   and  2,   respectively,   r-,  ^   is  the  distance  between  electrons 
1   and  2,   and  where   dx  ,  y   =  1,    2,   denotes  integration  over 

the  three  spacial  coordinates  of  electron  y.      These  Integrals 

LM 
will  be  expressed  in  terms  of  certain  auxiliary  functions   U 

LM 
and  V_  (a)   whose  calculation  is  described  in  the  appendices, 
P 

In  addition,  many  one-electron  integrals  can  easily  be  expressed 

LM        LM 
in  terms  of  U   and  V_  (oo  )   and,  at  the  end  of  this  section, 
a        p 

such  expressions  for  the  overlap,  kinetic  energy,  and  nuclear 
attraction  Integrals  will  be  given. 

We  first  note  that  only  the  factor  exp(im0)   depends  on  the 
sign  of  m   in  a  DO  and  that,  therefore,  a  DO  with  m  ^  0   is 
simply  the  complex  conjugate  of  a   DO  with  m  ^  0.   We  assume 
that  our  basic  set  consists  of  a  set  of  DO's 

X.  =  [2it)~   /^M.(^L)A.(A)exp(+im  0)  (3-2) 

with  m.  ^  0.   Thus,  a   j   denotes  a  pair  of  complex  conjugate 
DO'S   if  m.  >  0  and  a  single  DO   if  m  =  0. 

tJ 

The  volume  element  for  Integration  in  prolate  spheroidal 
coordinates  is 

dT  =  (R/2)^(A^-H^)dAd[xd0.  (3-3) 
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We   define 

O.  .(A,n)    -  M    (n)M    (n)A    (A)A.(A)(a2_p.2)  (3.4) 

and  may   then  write 

XX.dT    =    (2Tr)~^(R/2)^0.  .(A,^L)exp(lM0)di^dAd0  (3-5) 

to  represent  a  list  of  as  many  as  four  one-electron  charge  distri- 
butions differing  only  In  the  choice  of  signs  In 

M  =  +  "^-  ±  n^  •  •  (3-6) 

A  list  of  as  many  as  sixteen  J-lntegrals,  differing  only  In  the 
choice  of  signs  In  the  DO's  can  be  written 

J  =  Al^jl  /^l2^-2  ^i2^^ld'^2 

=  (2^r)"^(R/2)^/o^^.-^(A^,^l^)0^2j2^^2'^'2^'^''P^^'^l^l'^^^2^2^^ 


-1 
12 


X  dA^  dAg  dp.^  dM.2  d0^  d^^,  (3-7) 

where  M   and  M   are  the  M-values  given  by  (5-6)  for  the 
electron  1  and  2  distributions .   Noting  that   r-.^   depends  on 
0-,   and  0p  as  a  function  of  the  difference 

0  =  01  -  02.  (3.8) 

we  replace  0-,  by  0  as  a  variable  of  Integration  and  Integrate 
with  respect  to  0p.  The  result  Is  that  the  non-zero  J-lntegrals 
may  be  written, 

j(lll2,Ji,J2'^^"^^^^~^^^/^^y*^1^2^^P^^^^^^12^^1^^2^^1^^^2^'^'     ^^-^^ 
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where  M  =  M  =  -  M  .   (An  abbreviated  notation  Is  now  used  for 
the  Q's.)   Since   r-,  „   Is  Independent  of  the  sign  of  0, 


J(l-^l2J^J2.I^)  =  Ji^-^,^2'h'^2'~^^ 


(3.10) 


Therefore,  we  need  only  consider  M  >_  0  and  for  each  set  of 
Indices,   l-,ipj-,j',   there  are  at  most  two  Integrals  to  be 
calculated,  with  these  two  distinguished  by  the  choice  of  sign  In 


M  =  iTij.  +  m^  >  0 


(3.11) 


where  m^   and  m^   are  the  larger  and  smaller,  respectively, 

of  m.    and  m.  . 
il        Ji 

We  use  the  Neuman  expansion* 


-1 
\2 


2 
R 


1 


oo        M=L 
(2L+1) 


L=0 


I 


(-1) 


M 


M=-L 


(L-|M|)! 


( L+ I M I ) I 


pi^l(AjQl^l(A. 


.Pl^l(^,)pj^l(^2)e-^^ 


(3.12) 


M  M 

where  the   PJ  ' 's   and  QJ  ''s   are  the  associated  Legendre  functions 

of  the  first  and  second  kinds  and   A^   and   A^   are  the  smaller  and 

larger,  respectively,  of   A-,   and   Ap .   When  we  substitute  (3-12) 

in  (3.9)  and  integrate  over  0,   the  only  non-zero  terms  in  the 

summation  over  L  and  M  are  those  for  which  the  latter  is  the  same 

as  the  M   in  the   J-lntegral  (3.9).   Thus, 

-   °°         (L-M)!   ^ 


•  R, 


J{l^,l^,j^,3^,m)    =      (|)   ^  (2L+1) 


(L+M) ! 


■LM 


(3.13) 


See  reference  18,  page  l474  for  a  derivation  of  this  formula, 
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where 

"^  (L+M)!  "^l       1      ^  ^  112^12 

g^'^(A)    =  /      P^(^i)0    (A,n)dn,         7    =    1,2  (3.15) 

7  _-^         L  7 


The    Integral      T^^      is    transformed  by  a  method  due    to   RUdenberj 
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as  follows:   We  express   T^.^  as  the  sum  of  two  Integrals,  one  over 
the  region  where   A-,  ^  Ap  and  the  other  over  the  region  where 
A,  >•  Ap.   Then,  letting   A  =  A^    and   A '  =  A^   in  each  integral, 
the  result  is 


^LM  =  (-1^ 


(L+M)!  \  ^         {^^  ^  ^  ^  ) 


where 


G   (A)  =y    Pl(A  Jg   (A  JdA  .  (3-17) 

This  can  be  written 

T_  =  (_i)NiL,M)l  f    4^f  ^  Gf(A)G^^A)  /  dA .     (3-18) 
(L+M)!  "^1    P^(A)  (  dA    1      ^      ^ 


LM 


Integration  by  parts  yields 

.  (_,)M+1  (L,M)1  /-  U        ^Li^  /  af  (A)G^^(A).       (3-19) 
^^  (L+M)!  1    (dA   P^(A)  ^   ^      2 

The  formula 

P^A)  ^  Q?|(A)  -  Q^A)  ^  P^A)  =  ^-^^'''"'^^^^^•'  (3.20) 

^     dA   ^        ^     dA   ^       (A^-l)  (L-M)! 

gives 
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oo  G   (AJG   (A) 


Thus,  we  get 


oo       ,,  „,,  2°  g"(a)gM(,) 


J  =  (|)  V-  (2L+1)  -^ii=^  /   ^^^ ^ ^  dA.      (3.22) 

2   ^^j^         (L+M)!  1    (a2-1)P^(a)2 

From  the  definitions  of  G^^     in  (3-17),   g^^   in  (^.l^),  and 

n   in  (3-^),  we  obtain 

A  1 

Gy^(A)  =/  Pl(A')/\.^(A')Aj.^(A')  /  P^(^L)M.^(^x)]VI.^(^i)(A2-^x2)d^xdA' 
(3.23) 
Letting  , 

u^  =  /pf(n)n«M,  (^i)M,(^L)d^L,  (3.24) 


V^(a)  =/p^(A')A'\(A')A.(A-)dA',  (3.25) 

we  get 


p        -1 


g"^(>,)  =U«vLM(,)  _uLMvM(,).  (5.26) 

The  most  important  one-electron  Integrals  required  in  the 

calculation  of  the  wave  functions  and  energies  for  a  many-electron 
diatomic  molecule  are  the 

overlap,  * 

^  =   Ai  ^j^^  ^^-27) 
nuclear  attraction, 

\  =   -     S^l   -a"'  ^/-  ^5-28) 

\  =   -  /^i  ^b'   ^j^^  (5.29) 
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and  kinetic  energy, 

T  =  ft{-   |A);?f^.dT  (3.30) 

integrals.   From  the  definitions  given  above,  it  is  seen  that  for 
the  non-zero  Integrals  (3 -27 ) -(3 -30) j  one  must  have  M  =  0  and 
that 

S  =  (R/2)^G°°(oo)  (3.31) 

V^  =  -  (R/2)2[ug°V°°(oo)  -  U°°V°°(oo)]  (3.32) 

V^  =  -  (R/2)2[U°°V^°(ao)  +  uJ°V°°(cd)].  (3-33) 

Simple  formulas  for  the  kinetic  energy  integrals  in  terms  of  the 
nuclear  attraction  and  overlap  integrals  can  be  obtained  by 
multiplying  the  differential  equation  (2.1 )  for  X-       (o^  X  ■) 
by    X .      (or   X-)   and  integrating.   This  gives 

T  =  -  r  .V  -  Cv-V,  +  e.  s  (3.3^) 

^ai  a   ^bi  b    X  \^  >  / 

or 

T  -  -   ^a/a  -  ^bA  +  ^j  S-  (3.3V) 

These  two  formulas  with  i  7^  j   provide  an  excellent  check  on  the 
accuracy  of  the  one-electron  integrals. 

To  complete  the  discussion  of  the  selection  and  ordering  of 
the  non-zero,  non-identically-equal  two-electron  integrals,  it 
remains  only  to  describe  the  selection  and  ordering  of  the 
i-i  j'-i  ip  Jp  's-   First,  we  form  a  list  of  distribution  functions 
by  listing,  in  numerical  order,  all  pairs   i,  j   with   i  ^  j,   and, 
with  each   i,  j,   we  list  the  one  or  two  possible  values  of 
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M  =  m,^  +  m^ .   Each  triple   Cl,j,M)   then  corresponds  to  a  set 
G^^(a),  L  =  M,  M  +  1,,..   of  G-functions  while  the  triples 
with  M  =  0   correspond  to  the  one-electron  integrals  described 
above.   Then,  for  each  distinct*  pair  of  triples   (i,  ,j-,,M,  ) 
CipjJp^Mp)   having  M-,  =  Mp  we  use  the  corresponding  pair  of 
sets  of  G's   in  formula  (3-22).   In  the  case  of  the  homopolar 
molecule,  there  is  one  more  selection  rule:  the  one-electron 
integrals  above  are  zero  unless   ^.p  +  ^.p   is  even  and  the  two- 
electron  integrals  are  zero  unless  H .  ^    +   jl  .-,      is  of  the  same 
parity  as  i^^   +  i .p. 

The  above  procedure,  used  by  Harris   and  others,  enables  one 
to  store  and  save  G-functions  for  all  pairs  of  orbitals   /(  •   'X  • 
and,  during  the  variation  of  non-linear  parameters,  alter  only 
G-functions  involving  orbitals  whose  parameters  have  been  altered. 
This  represents  a  considerable  economy  in  computation  as  compared 
with  other  methods  requiring  a  complete  calculation  of  all  inter- 
mediate functions  in  each  two-electron  integral  for  which  at  least 
one  of  the  four  orbitals  has  been  altered. 
4.    Concluding  Remarks 

The  ultimate  evaluation  of  the  usefulness  of  DO's   as 
compared  with  STAO's   or  simple  spheroidal  coordinate  orbitals 
can  be  made  only  after  evaluating  the  accuracy  of  computed  results 
in  terms  of  the  computational  effort  needed  to  obtain  them.   In 
the  present  formulation,  the  machine  program  for  the  IBM  7090 
takes  only  40-50  seconds  to  calculate  all  one-  and  two-electron 

Distinct  pairs  of  triples  are  obtained  and  ordered  by  requiring 
that  either   i-,  >    ip   or  i-,  =  ip  and   j'  >   j'  . 
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Integrals  for  the  DO's   Is,   2s,   3s,   Jdcr,  2pcy,   ^po^,   ^po',  hfv, 
2p7r,   and  3d7r,   (giving  485  two-electron  Integrals)  for  parameters 
appropriate  to  the  Hp  molecule  at  R  =  1.4,   the  ground  state 
equilibrium  distance.   This  is  less  time  than  is  required  for  some 
of  the  STAG  integral  programs.   However,  this  is  more  time  than 
would  be  required  for  the  spheroidal  coordinate  orbitals   Xpo 
since  the  same  methods  are  used  for  both  and  the  latter  are  simpler 
functions.   An  important  incentive  to  the  use  of  DO's   would  result 
if  the  computational  effort  for  integrals  of  DO's   could  be  reduced 
by  improved  integration  methods.   Perhaps  one  could  find  some  method 
taking  advantage  of  the  fact  that  the   DO's   and  their  factors 

satisfy  simple  second-order  differential  equations.   Some  of  the 

17 
ideas  in  the  method  used  recently  by  Auffray,  Percus,  and  Ross 

for  atomic  orbital  integrals,  in  which  the  two-electron  integrals 

were  transformed  to  integrals  over  momentum  space,  may  be  applicable 

here . 
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Appendices 

I.   Computation  of  the   DO's 

In  the  present  calculations,  only  the  homopolar  molecule,  with 
P-,  =  0,   is  considered  and  the  solution  of  (2.'^)  is  expressed 


oo 


k=mT6  ^ 

Where   5  =  0(=l)   if   ^  +  m   is  even  (odd).   The  summation  is  in  steps 
of   2   over  all  even  (odd)   k,   with  k  ^  m,   if  i      is  even  (odd). 
Substituting  (l.l)  in  (2.4)  gives  the  recurrence  relations 

p   =  F   =0 
-2  -1     ' 

A,  F,-    +  B.  P,  +  C,  P,  ,^  =  0  (1.2) 

k  k-2    k  k    k  k+2 

where , 

(k-m) (k+m) 


\  = 


(2k-l)(2k-5) 


B   =  -k(k+l)-A  _^    (k+m)  (k-m)  ^    (k+m+l)  (k-m+l) 
^^       p^     (2k+l)(2k-l)     (2k+l)(2k+3) 

^   _  (k+m+2) (k+m+l) 
^k  -  • 


(2k+3)(2k+5) 

The  eigenvalues  are  those  discrete  values  of  A  for  which  the 
series  in  (l.l)  converges  for   -  1  Jl  M-  ^  1  with  Fj^'s   satisfying 
(1.2) .   This  gives 

W2  ;^(p/2k)^  (1.3) 

for  large   k. 
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The  form  for  the  solution  of  (2.5)  used  here  Is 


m 
A(a)  =  (A^-l)^e"P^(A+lfy^  G^^z'^  (1.4) 

k-0 

where  C  =  PoP"   -  m  -  1, 

z  =  (A-1)/(A+1) 
and  the  ^  's  satisfy 

A,  G,  .  +  B,  G,  +  C,  G,  ^,  =  0,  (I.^) 

k  k-1    k  k    k  k+1     '  ^      ^' 

with 

A^  =  (k+m-p2/p)Ck-p2/p) 

B^  =  A-p^-2pm-2k^-(l+2k) (2p+m+l)+(p  /p) (2k+2p+m+l) 

C^  =  (k+l)(k+m+l) . 

Convergence  of  the  series  In  (1.4)  with  the  recurrence  relations 
(1.5)  gives 

G   /G   ->  1  -/is 

'^k+i/^k  --  '  y  1, 

for  large   k. 

Both  recursion  relations  (I.2)  and  (I.5)  a^e  essentially  of 
the  form 

«k  Vi  +  Pk"k  +  ^k^k+i  =  0  ^^-6) 

with  the  asymptotic  behavior  of  the    ^v-'^'  ^°^   large   k, 
expressible  In  the  form 
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^k  =  ^k+1  \  (1-7) 

where  R,   is  known.   To  make  the  present  method  applicable  to 
a  more  general  class  of  three-temn  recursion  relations,  we  will, 
for  the  moment,  merely  assume  that  the   a,  ,  p,  ,   and  7i^'s   and 
their  derivatives  with  respect  to  the  eigenvalue  parameter  A  are 
computable.   The  procedure  is  to  assume  a  value  of   A,   give  a  small 
arbitrary  value  to  ^m. -i  >      where   N  is  a  sufficiently  large  integer, 
use  (1.7)  to  obtain   x   and  apply  (I.6)  for  k  =  N,  N  -  1, . . . ,  m, 
where   k  =  m   is  the  first  index  for  which 

la  X   -,  I  ■*=  I7  X  ,-,  i  .  (1.8) 

I  m  m-1 1  —  r  m  m+1  I  ^ 

The  inequality  (I.8)  implies  that  greater  numerical  accuracy  can 

be  obtained  for   x  ■^  >      x   p^  •  •  •   by  recurring  outward.   An  arbitrary 

value  is  then  given  to  Xj-,  and   (I.6)   is  used  with  increasing  k 

up  to   k  =  m.   Letting  the  x,  -values  of  the  latter  step  be  denoted 

by   Yq,    y-^,...,  y^^  and  letting 

1^   -  Xj^y^^  x^      for  k  =  m  +  1,  m  +  2,.  ..,  N    (I.9) 

we  thus  obtain  a  trial  solution  yQ,  y.,  ,  .  .  . ,  yvr, -i   satisfying 
(1.6)  for  all   k  ;^  N  except   k  =  m.   The  amount  by  which  the  m-th 
equation  of  (I.6)  is  not  satisfied  is  regarded  as  a  function  of  the 
trial  eigenvalue  A, 

F(a)  =ay  -,+py  +7y,-,.  (1. 10) 

m'^m-l    ^m-'m   ^m-'m+l 

The  Newton-Raphson  method  is  used  to  find  the  zeros  of  F(A) .   By 

this  method,  successive  Iterates,   A  ,  A  ,...,   converging  to  a  zero 

of  P(A)   are  given  by  the  formula 

A^+1  =  A^  -  F(A^)/F'(a'').  (1. 11) 
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Letting  primes  denote  derivatives  with  respect  to  A,   we  can 
differentiate   (I. 10)   to  get 


F'  (A)  =  a'y  .  +  B'y  +  7'y  -, 

m-^m-1    ^m*^m   ^m"^m+l 


+  ay'-,+6y'+7y'.  (1.12) 

The   Yi^'s   are  defined  as  dlf f erentlable  functions  of  A  by  the 

procedure  described  above.   The  derivatives  x'   ^,   x'    and 

m+1 '    m 

y'   are  calculated  by  using  the  recurrence  relation  obtained  by 
differentiating  (1.6)  v;lth  respect  to  A.   With  starting  values 
y'   =  x«  =  0,   the  calculation  of  the  ^A's   and  yA's   is 
performed  simultaneously  with  that  of  the  ^i^'s   and  yi^'s.  Finally, 
y  •  -,   Is  obtained  from  the  formula  obtained  by  differentiating  (1.9) 
with  k  =  m  +  1,   giving 

y',^  =  (x'  ,y   +  X   ,y'  -  y   ,  x')  x""'"  .  (1.13) 

''m+l     m+l*^m    m+l-^m   *^m+l   mm  v   -^v/ 

The  above  procedure  yields  a  finite  series  approximating  the 

infinite  series  for  the  exact  solution.   It  has  been  shown   that, 

for  the  equations  above,  the  successive  ratio  of  any  two  expansion 

coefficients  can  be  represented  as  an  infinite  convergent  continued 

fraction  and  that  the  condition  for  an  eigenvalue  is  that  the  continued 

fraction  be  equal  to  the  same  ratio  as  calculated  by  recurring  outward 

from  the  first  term.   Since  the  procedure  described  above  yields 

expansion  coefficients  whose  successive  ratios  are  the  same  as  would 

be  given  by  the  continued  fraction  procedure,  we  can  conclude  that 

if  we  take  N  sufficiently  large  we  can  approximate  the  eigenvalue 

and  any  finite  set  of  expansion  coefficients  to  any  desired  degree 

of  accuracy. 
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Subroutines  for  solving  for  A   and  A,   and  the  expansion 
coefficients  of  the  corresponding  eigenfunctions   are,  in  the  machine 
program,  used  in  an  iterative  procedure  which,  by  successive  interpo- 
lations, finds  the  desired  zero  of 

A  (p)  -  A^(p)  (1.14) 

for  each  DO  given  to  the  program. 

The  expansion  coefficients  and  other  wave  function  parameters 

of  the   M(ij.)   and  A(a)   solutions  used  here  have  been  tabulated 

"1  ft 
by  Bates,  Ledsham,  and  Stewart   for  a  number  of   DO's.   These  were 

Invaluable  in  checking  the  machine  program  and  for  providing  first 

approximations  of  p  and  A   in  the  iterative  procedures  described 

above.   The  present  program  provides  considerably  more  accuracy  than 

is  given  by  Bates  et  al .  and,  in  the  numerous  cases  tested,  was 

found  to  agree  with  the  tables  to  the  degree  of  accuracy  claimed  by 

the  authors . 

LM 

II,   Calculation  of  the  U    -  Functions. 

a 

The  following  describes  a  method  for  calculating 
1 
U^^  =  /  Pf(n)n>.(n)M.(^x)dn  (II. 1) 

where  M.(|j.)   and  M.Cia.)   are  of  the  form  (l.l).   Let 

6.  =  0(l)   if  I.    +   m.   is  even  (odd), 
1  11 

and  (II. 2) 

A  =  0(l)   if  L  +  M   is  even  (odd). 

Here,  and  in  what  follows,   i   and   j   indices  will  appear  with 
the  wave  function  parameters  of   X.   and  X.,,      respectively. 

J.  J 
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It  is  assumed  that   i   and   j   are  assigned  so   that  m.  ^m.. 

We  define  the  Integrals 
1 

FJ(a)  =  /  Pf-iu-)ix^n^ii^)d^.  (II. 3) 

Then, 

,        .   „    (k+m. ) ! 

F,^(0)  =  Fi  -^ i_   .  (II. 4) 

^       ^  2k+l  (k-m^) ! 

where  the  F?;'s   are  the  expansion  coefficients  of   (l.l)   for 
M.(ia).   From  (2.4),  with   p-j^  =  0,  we  get 

H^M.Cn)  =  p-2[j:^^(^L)+A^]M.(n)  (II. 5) 

which,  when  substituted  in  (II.3),  gives 

FJ(2)  =  p-2[k(k+l)+A.]FJ(0)  (II. 6) 

Next,  we  define  F^(K',a)  by 

^  /  .      m . 

(l-n2)^/2p.«M.(M.)  =  Yi        Pi^K,a)P^^{ii)  (II. 7) 


k 
where 


K  =  m.  -  m_.  >  0  (II. 8) 

and  the  summation  is  in  steps  of   2   over  even  (odd)   k's   if 
a  +  i.  +  m.  +  m.   is  even  (odd).   Then,  obviously, 

cj         -^         J 

pJ(o,o)  =  P,-^"  (II. 9) 

k   '      k 

and,  by  the  method  used  to  derive  (II.6),  one  gets 

FJ(0,2]  =  p:^[k(k+l)  +  A  ]fJ(0,0]  (      ;     (11.10) 
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From  relations  satisfied  by  the  Legendre  functions,  one  gets 
recursion  relations  for  raising  the  Indices  K     and  a: 


-1t,J 


-1, 


P^{K,a)    =    (2k-l)   Fj[_i(A-l,a)  -  (2k+3)   Fk+l^*^"^'"^       (II. ll) 

k-m.   .             k+m.+l 
Fi{K,a)    =    ( -)'Pi   AK, a-1)    +    ( i — )F,^_^AK,a-l)  (11.12) 


k' 


2k-l 


k-1 


2k+3 


k+1 


The  F>'(/c,a)   with   a  ^  2  are  obtained  by  using   (II.  // )   and 
(II.  I2)      with  the  starting  functions   (II.'^}   and  those  with  a  *=  2 
are  obtained  from  the  starting  functions  (II. 8). 
Now,  we  write 


uKX  _  (2K)! 
"   "2V. 


f    M.(m,)[(1-h2)^/2^^j(h) 


d|a 


rn"Mi(n)[(l-n^)^/^^L"~^Mj(n)' 


If   a  <  2 


d\i        If   a  >  2 


(11.13) 


Substituting  the  series  (II. 7)  for  the  quantity  In  brackets  gives 


U 


/      \           Fi(0)PJ(X,a) 
/ok            k 

a. 

_    (2X)! 
2^K! 

•   X 

,  k=m . +5 . 

/      >     2      F^(2)PiJ('<,c-2) 

.  k=m . +5 . 
\        11 

If   a  <  2 


(II. 14) 


If   a  >  2 


One  may  write  M   In  the  form 

M  =  2t  +  K 
where   t  =  0  and  m.   for  the  two  positive  values  of  M  In  the 

J 

choices  of  signs  In  (3.6).   For  t  =  m .  f^  0,   we  obtain  U    by 
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m .   applications  of  the  recursion  relation 

uM+2,M+2  ^  (2M+3)(2M-M)(U^'^-U^,).  (ll.l-) 

a  a   a+d 

LM 
Finally,  we  obtain  U    from 

U^^  =  (L-M)"^[(2L-1)U^:|'^-(L+M-1)U^-^'^].      (II.16) 
a  a+l  a 

LM 
The  non-zero  U   's   are  those  for  which  a  +  L  is  of  the  same 
a 

parity  as  £.+£..      For  the  two-electron  integrals,  we  only  need. 

LM 

U   with   a  =  0  and  2.   Therefore,  the  sum  in  (5-22)  is  only 

over  those   L's   having  the  same  parity  as  £.+£.. 

III.  Calculation  of  the  V   (a ) -functions 

The  calculation  of  the  functions 

A 
V^(A)  -/PL(A')A'\(A')Aj.(A')dA'  (lll.l) 

is  performed  by  using  the  recursion  relation  (II.16),  with  U's   and 

a's   replaced  by  V's   and   P's,   and  starting  with  the  functions 

A       M 

vf(A)  =^/(A'2_i)2,.Pa  (A')A,(A')dA'.    (ill. 2) 
P        2  1^!   1  "^ 

We  make  the  change  of  variables 

z  =  (A-1)/(a+1)   ,   z'  =  (a'-1)/(a+1)  (III. 3) 

LM 
and,  letting  Vg   be  regarded  as  a  function  of  z,   we  have 


0 
where 


vT'(z)  =   /f(z<)dz'  (III. 4) 


^MM, 


M 

^•(2')  = -4|t^  (a'2-i)2a.Pa.(a')A,(a')(a'+i)2  (III. 5) 

2    M! 
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Letting 

h   =  1/N 


^k  =  ^^ 


A^  =  (l+Zj^)/(l-Zj^)  (III. 6) 

^k  =  ^^^k^ 

dz 

MM 

we  get  Vg  (z,  )   at  the  equally-spaced  points   z,   by  use  of  the 

Integration  formula 

The  derivatives  in  (ill. 7)  are  calculated  with  the  formulas  obtained 
by  differentiating   f(z'). 

The  above  method  yields   6  and  7   significant  figures  in  the 
two-electron  integrals  with  50  integration  points  for  integrals 
appropriate  to  the   H_  molecule  at  the  internuclear  distance 
R  =  1.4  a.u.   Improvements  being  considered  are  to  replace  (ill. 5) 
by  the  change  of  variable   z  =  1/a   and  to  use.  Instead  of  (ill. 7) 
the  higher  order  integration  formula 

^p  ^^k+1^  =  ^p  ^^k^  +  a^V^k+i^  +  To^^k-^k+1^ 

^J2^(fn+fH   )  (III. 8) 


120^"k   k+1 

Simple  expressions  for  the   f'v's   can  be  obtained  by  using  the  fact 
that  the  f\i'K)      factors  in   f(z)   satisfy  the  second  order  differ- 
ential equation  (2.5)- 
IV.   Calculation  of  T^^j^ 

The  above  procedure  yields   G   (a)   at  equally-spaced  z-points 
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The  integral 


T 


LM 


.°°  Gf(A)Gf(A) 
1    (a2-1)P^(a)2 


dA 


(IV. 1) 


is  transformed  to  an  integral  over  z, 

1 


^LM  =^F(-)dz, 


( IV . 2 ) 


where 


F(z)  = 


Gf(A)Gf(A: 
2z  P^(A)2 


(IV. 3 


and  is  integrated  by  the  same  formula  as  was  used  in  (ill. 7)-   The 
resulting  integration  formula  for  T^^^   is  however  considerably 
simpler  since  the  derivative  terms  cancel  at  all  interior  points 
giving  the  formula 


T   =  h 
LM   " 


iF(o)  +  ipd)   + 


N-1 


,  2  - 

n 


2—     ■    2--'  ■  l_^^^k^J  ■  12 


k=l 


F' (O)-P'(l) 

(IV. 4) 


where  primes  denote  derivatives  with  respect  to  z .    To  evaluate 
the  integrand  and  its  derivative  at   z  =  0,   we  first  define 


A(z)  =  (A+l) 


p/p-m^-pA 


OD 


F"  V 


k 


k=0 


so  that 


Then, 


m 
2  .  n2, 


-1" 


A(A)  =  (A^-l)^(A+l)--'A(z) 


(IV. 5) 


(IV. 6) 


29 


m.+m  . 

2 i 


G^(A)  =  |/p^(a)(Uqa2-U2)(a2-1)   2   /^  (^  )7^^.  (^  )dz 

(IV. 7) 


m.+m. 
^m. +m  .   ^_^  ^  ^  ^ 

z 


m.+m  .+2 

1   J 


~  M  ~I 

Pt-(i)W.  .+  higher  order  terms  in  z 

_  L'    ij    ^  J 

(IV. 8) 
where 

^ij  =    (Uo-Ug)^  Ai(0)A.(0),  (IV. 9) 

the  notation  denoting  that  U^-Up   Is  calculated  with  M.(^l) 

and  M  .(|j.)  . 

Then, 

22e-l^e+l 

^^^^  =(m.,+m.,+2J(m.^+m.^+2)^^lljA2j2  +  higher  order  terms  In  z] 
il   jl      i2   j2        '^ 

(IV. 10) 
where   e  =  ^{m^-^+m.-^+m^^+m.^) 

Therefore, 

F(0)  =  0 
and  A 

^^°^  -  V  (IV. 11) 

V  0  otherwise . 

It  can  easily  be  shown  that,  at  z  =  1,   P(z)   and  Its  derivative 
are  given  by 

F(l)  =-F'(l)  =  I  G°°(oo)G°°(oo)  (IV. 12) 

* 

If   L  =  0  and  are  zero  otherwise.   Then,  we  write  the  Integration 
formula  as  the  simple  vector  product 
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T,..  =V^gH'  G^:'  tlV.13) 


1 


where 


LM  /_ Ik  2k 

k=0 


,LM 


GrCAkI 


G^M^/h     7 ^     k=l,2,...,N-l,  (IV. 14) 


l^"/2-k    P^(Ak) 


G^ 
7O 


-^  W.   .    if  m..   +  m..^  ?^  0  (IV. 15) 

i^/g   17,  J7       17     37 


=   0  otherwise 


|x/h  G°°(qd  )    If   L  =  0 
g"^    )        ^  (IV. 16) 


7N 


=  0  otherwise 


.00/ 


h/2y^G;^(co)   If  L  =  0 

G^M   ^  )  (IV. 17) 

7,N+1  ) 

=  0  otherwise 

In  the  machine  program,  the  quantities   G^,  ,  L  =  M  +  A,  M  +  A  +  2, 

and  k  =  0,1,...,N+1   are  calculated  and  stored  on  magnetic  tape  for 

V  =  1,2,5,...,   where   v   is  an  Index  going  over  all  pairs  of  DO's 

X.H.,    1  ^  1.   The   v's   are  partitioned  into  sets  with  the   v's 
1  J    - 

of  a  set  characterized  by  having  the  same  value  of  M  and  A. 
The  program  then  goes  over  all  pairs,   v, ,  Vp  with   v^  >  v^   in 
each  set,  forms  the  vector  products   (IV. 15)   with  L  =  M  +  A  , 
M  +  A  +  2,  ...    and  uses  these  in  formula   (3-13)   to  yield  a 
J-integral.   Prom  the  discussion  at  the  end  of  section  5  above,  it 
is  seen  that  this  procedure  yields  all  non-zero  two-electron  Coulomb 
repulsion  integrals . 
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